Critical behavior of the Ising model with long range interactions 
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We present results of a Monte Carlo study for the ferromagnetic Ising model with long range 
interactions in two dimensions. This model has been simulated for a large range of interaction 
parameter a and for large sizes. We observe that the results close to the change of regime from 
intermediate to short range do not agree with the renormalization group predictions. 
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While the Ising model long range interactions (LRI) 
has been studied for a long time as a generalization of 
the Ising model with short range interactions, very few 
results have been obtained for the case of weak inter- 
actions decaying faster than the dimension d since the 
first works in the 70's We will recall some of these 
results after defining the model that we will consider in 
this letter. The Ising model with LRI is defined by the 
Hamiltonian: 



-R = 



J 



<%3> *J 



(1) 



with spins taking values S'i = ±1 on the sites i of a regular 
lattice and Vij the distance between the spins on the sites 
i and j . The sum < ij > is over all the pair of spins and 
we will consider only a ferromagnetic interaction J > 0. 

In 1972, Fisher and al. [2^] performed a renormaliza- 
tion group (RG) study for the 0{n) model with LRI (the 
Ising model corresponding to n = 1). They identified 
three regimes (for d > 1) : i) the classical regime with 
a < d/2 which is believed to be with a mean- field behav- 
ior, ii) an intermediate regime for d/2 < u < 2, iii) the 
short range regime for a > 2. The value a = d/2 marks 
the border of the mean field regime, meaning that the 
ordinary perturbation parameter e = 4 — d is replaced by 
2(7 — d. In the relation rj = 2 — a was also conjectured 
in the intermediate regime. This result was questioned 
since for tr = 2, the exponent r] vanish while for a > 2 
its value has to be the one of the short range model r]sr 
which is positive for d < 4. Then it would imply a jump 
of the exponent 77 from up to rj^r at a = 2. This point 
was first considered by Sak [3|, who, by taking in account 
higher order terms in the RG calculations, predicted that 
the change of behavior from the intermediate to the short 
range regime takes place at cr = 2 — rjsr- Many other 
studies have considered also this problem with various 
conclusions. In particular, van Enter [3l obtained that 
for n > 2, long range perturbations are relevant in the 
regime 2 — rjsr < ct < 2 in contradiction with Sak results. 
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Gusmao and Theumann J^, by considering a develop- 
ment in terms of e' = 2a — d in place of e = 2 — cr, 
obtained a similar result, namely the stability of the long 
range perturbation for a < 2. 

Note that all these studies are using a renormalisa- 
tion group approach with an e expansion of a Landau- 
Ginzburg effective Hamiltonian such that the propagator 
contains a p'^ term in addition to the ordinary term. 
While it is know that this approach gives accurate results 
for two and three dimensions for the short range model 
(with just the term) 0, it is only after comparing 
these predictions with other methods, numerical, high 
temperature expansions, (and the exact result in two di- 
mensions), etc. that we can believe in these predictions. 
Similar comparisons need also to be done when consid- 
ering the case with LRI and this is the main purpose of 
the study presented in this letter. 

A first numerical study of the exponent rj for d = 2 
as a function of a has already been done by Luijtcn 
and Blote In particular, they obtained in the inter- 
mediate regime a result well described by the exponent 
rj = 2 — cr up to 2 — (T = Tjsr and rj = rjsr = 1/4 for 
larger cr. Thus their measured exponent is in agreement 
with 7] = max(2 — cr, 1/4) which corresponds to the RG 
predictions of [3|. 

In the present study, we improve Luijten and Blote 
study. In particular, we repeat the measurement of 77 
close to the region where its behavior is changing, i.e. 
for cr ~ 2 — rjsr ■ We confirm their result that there is no 
discontinuity but we measure a clear deviation from the 
behavior predicted by Sak 

In order to be able to consider large lattices, we need 
to employ an efficient algorithm. Since the model is fer- 
romagnetic, we can employ a cluster algorithm which will 
reduce the auto correlation time, that is the number of 
upgrades between two successive independent configura- 
tions. For the short range Ising model, Wolff algorithm, 
which builds a single cluster per update, is the most ef- 
ficient one We will adapt this algorithm to the case 
with LRI. The algorithm can be summarized by the fol- 
lowing steps. We start from a spin Si at some randomly 
selected position i. Next, we add to this spin any other 
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spin Sj with a probabihty 



the process until we obtain 
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with /? the inverse temperature. We repeat the same 
operation with all the added spins. In order to build a 
cluster containing M spins on a lattice of N sites, we have 
to compute the probability 1^ ior ^ N x M bonds, then 
the number of operations is 0{N x M). This number of 
operations will be drastically reduced by our algorithm 
that we will describe now. We will present here only the 
main steps of this algorithm, a more complete version 
will be presented elsewhere J9|]. 

The main idea is that since p{rij ) is very small for large 
rij, then it is much faster to compute the probability 
that one (or more than one spin) among all the spins 
at this distance are connected to the original spin Si. 
In practice, to build this cluster we proceed as follow. 
First, starting from some arbitrary position i, we order 
in some pile Vi{k = 1, • • • , ~ 1) of length — 1 all 
the other positions on the lattice, ordered in function of 
the distance. To be more precise, the pile will contain the 
difference of positions between the point i and j . The pile 
■pi will be the same for any i, so this operation is done 
only once. Next we consider all n{r) spins at a distance 
r. This is a very fast operation if we build a second pile 
7^2 (^) containing the position of the first spin in Vi at 
a distance r. Then n(r) = V2{r') — 'P2{r), with r' the 
smallest distance on the lattice such that r' > r. For a 
given configuration of spins and a given starting point i, 
we need only to consider the spins with the same value 
as Si and which are not already part of the cluster under 
construction. We denote by n'{r) < n(r) the number of 
such spins. Then the decomposition 



1 = 



((l-p(r))+p(r))"'W (3) 
p(r)"'('') + n'(r)p(r)"'('')-\l - p(r)) 
+ --- + (l-p(r))"'M 



fe=0 



will contain the probability of having zero spins 
connected p(r)" of having one spin connected 
n'{r)p{r)"' — p(r)), two spins connected 

n'lr){n'{r) - l)p(r)"' ^-^(i _ p{r)f /2, etc. Next 
we randomly generate a number e in the interval [0, 1]. 
If e < p(r)" (''^ no spin will be reversed at the distance 
r. We can then ignore all the n(r) spins at distance r, at 
the cost of generating one random number. Otherwise, 
if p{r)"'W < e < p(r)"'W + n'(r)p(r)"'W~i(l -PW), 
one spin has to be reverse among the possible 
7i'(r) spins. We then need a second random num- 
ber to choose one spin among the n'{r). Next if 
p{r^"- ('■) -|- n'(r)p{r)^ _ p{r)) < e, we continue 



(n'(r)~fc)!fc! 



with kfnax the number of spins that we have to reverse 
and thus the number of random numbers that we need 
to generate to choose these spins among the n'{r). 

For each distance r one needs also to compute n'[r). 
This corresponds to the delta function in eq.©. In 
principle, this means that we need to perform an ad- 
ditional number of operations n{r) to determine n'{r). 
In fact, in most cases, this part can be skipped. In- 
deed, it is much more convenient to first compare e with 
p{j-'^n(j-) ^ p^T-)" (»•), Thus we will compute n'{r) only if 
p(j,-jn(r) ^ which will be very rare for r large. Then in 
most cases one can consider all the n{r) spins in the slice 
at distance r with the cost of generating a single random 
number. 

A second improvement is to consider an ensemble of 
successive slices with ri < r < r2. The decomposition 
^ is then replaced by 



1 = X{{{l-p{r))+p{r)r'^^ 

r—7'i 



(5) 



p{r) 



The choices of ri and r2 are guided by efficiency. Start- 
ing from a first slice ri , we will add slices up to r2 with 
the condition that YirLr p('')"''^'' remains close enough 
to 1 to ensure that in most cases the randomly gener- 
ated number e will be such that e < nr=ri p('')"*-'^'' < 
JJr=ri Pi''')^ '"^^ ■ Again, we will first compare e with 
rirLri Pir)"^^^ ■ If e is smaller than this product, we can 
then skip all the spins in the slices of distance ri < r < r2 
for the cost of generating a single random number. Oth- 
erwise, we have to count the number of spins n'{r) in 
each slice ri < r < r2 which can be connected to the 
cluster {i.e. spins with the good sign and not yet part 
of the cluster). This corresponds to A^n.ra = X^ri 
operations. Now, if the condition nr=ri p(^)"^'^^ < e < 
nrLriP(^)"'^"^ is satisfied, we can skip all the spins in 
the slices of distance r with ri < r < r2 and with a cost of 
generating a single random number + Nr-^^r2 operations. 
Otherwise, we have to reverse at least one spin and we 
have still some additional operations. It can be checked 
that the number of additional operations in that case is 
smaller than AVi,r2- 

The choice of ri and r2 can now be obtained by re- 
quiring that the probability of reversing one spin times 
the number of operations is comparable to the number 
of operations for having to reverse no spins, i.e. when 
eq.Q is satisfied with kmax = 0. As said before, for 
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this case, the number of operations is C(l). As a first 
approximation, the probability of reversing one spin is 
P,.^ ,-2 ~ J2r=ri ^ p{^))- Then the optimal choice 

ri and r2 is such that Pri,r2 ^ ^ri,r2 — C'(l)- A simple 
extension of this argument Q leads to the result that 
the total number of operations for a system with N spins 
grows like 0{N). 

Our algorithm is similar in spirit with the one devel- 
oped by Luijten and Blote but the implementation is 
rather different. In their algorithm, the second part of 
the weight ^ is obtained by building a cumulative bond 
probability. In order to perform this step efficiently, they 
need to approximate the interaction between two spins 
by a integral. While this approximation does not af- 
fect the universal critical properties, other nonuniversal 
quantities like the critical temperatures will be different 
from the one that we obtained. It was argued in 
that Luijten and Blote algorithm requires 0{N log N)) 
operations. 

Note also that while we present here results for the 
case of the Ising model in two dimensions, our algorithm 
is valid for any Potts model and in any dimensions. We 
will present elsewhere results for the case of the Ising 
model with long range interactions in three dimensions 
i- 

We will now report on simulation done on a triangular 
lattice of linear size up to L = 5120 with periodic bound- 
ary conditions. To implement these boundary conditions 
with long range interactions, we employed the minimum 
image convention. In order to be able to obtain a good 
precision, we had to accumulate a lot of statisics. For 
each value of a,K = (3 J and size L, we accumulated 
statistics over 10^ x t updates, with r the autocorrela- 
tion time. The time for the update of one sample of size 
5120 X 5120 is ~ 0.1 second for a = 0.8 and ~ 1 second for 
CT = 1.8 on a recent PC. With the help of the Wolff clus- 
ter algorithm, the autocorrelation times are rather small 
even for the largest sizes considered. For L = 5120, we 
determined r ~ 200 for a = 0.8 and r ~ 66 for cr = 1.8. 
The total computing time required to produce the data 
presented in this work corresponds to 0(100) years of 
CPU time on a single core processor. In order to test the 
universality of our results and to have a more direct com- 
parison with the results of Luijten and Blote [3] , we also 
performed a simulation for cr = 1.75 on a square lattice. 

For each value of a, we first had to determine the crit- 
ical value K,.. This was done by considering a magnetic 
cumulant similar to the Binder cumulant ll| and which 
is an adimensional quantity. We will consider the mag- 
netic cumulant defined by B{L,K) : 



B{L,K) 



2\2 



1^ ' 



(6) 



with (m*) the thermal average of the magnetization to 
the power i. This cumulant will converge to the value 
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FIG. 1: (Color online) Exponent r] vs. o computed with data 
up to size L = 1280, 2560 and 5120. The inset contains a 
magnified part for 1.5 < cr < 2.2 with in addition the data ob- 
tained by Luijten and Blote 0. The dotted lines correspond 
to the prediction from the RG analysis, rj — max{2 — a, 1/4). 
The dashed line connecting our measured points is a guide to 
the eye. 



1 in the ferromagnetic phase and 1/3 for the paramag- 
netic phase. The curves describing the cumulant versus 
K for different sizes L will cross at a value Kc{L) which 
will converge toward the real critical point in the large 
size limit. At the crossing point, the cumulant B{L) will 
converge towards a finite value. For the critical Ising 
model with short range interaction, the limiting value 
was computed on the triangular lattice with a value of 
ymYL^^B{L,Kc) = 0.85872528(3) Close to the 

critical point, the finite-size scaling behavior is expected 
to be of the form 



BiL,K):^f{{K-K,)L'/n 



(7) 



with f{x) a dimensionless function. Then the crossing of 
the curve B{L, K) as a function of K for two different 
sizes L and L' takes place at if = Kc- In the following, 
we will always choose L' — 2L and express the computed 
quantities as a function of L only. In practice, due to 
corrections to scaling, we need to take into account addi- 
tional correction terms to B{L,K). We will consider in 
the following the leading correction of the form : 



B{L,K)^f{{K~K,)L 



(8) 



With such a term, the crossing of the curves will take 
place at a size dependent Kc{L) '■ 



K,{L) ^Kc + aiL- 



(9) 



with ai a constant and wi = yi + l/v. We will then de- 
termine the value of Kc by extrapolating the measured 
values of the crossing Kc{L) while including one single 
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correction. We have checked that with such a procedure, 
the values of wi and Kc are stable if one keeps only large 
size data, L > 100. The same will be true for the com- 
puted values of 77, see below. We have also tried to in- 
clude further correction terms in ([5]). We have checked 
that the addition of subdominant corrections does not 
affect the results for L^ax > 1280, i.e. the change in Kc 
is much smaller than the error bars on this value. We 
also performed a fit of B{L, K) close to Kc with a devel- 
opment in powers of [K — Kc)L^^'^ and with terms L~^' 
corresponding to correction to scaling. The obtained ex- 
ponents are always in good agreement with the ones from 
a direct fit of the form ([5]). In particular, we obtain that 
ly remains very close to 1. We quote 1^ ~ 0.96(2) for 
(T = 1.6. Finally, in 7], it was claimed that logarithmic 
corrections to scaling were needed in order to obtain a 
consistent analysis for a = 1.75. This is not the case 
in the analysis of our data. In fact, even if we impose 
the existence of logarithmic corrections, we observe that 
most of our results are not affected. These corrections 
will mostly affect the exponents yi without a measurable 
change on the exponent rj for the largest sizes that we 
can simulate Q. 

In the following, we will be interested in the region 
cr > 1.4 since this is where we obtain new results. Our 
main results are shown in Fig. [TJ It contains our nu- 
merical results for the exponent i] = 2/3/1/ versus a ob- 
tained for three ranges of linear sizes, L^ax — 1280, 2560 
or 5120. In each case, L^ax corresponds to the maxi- 
mum sizes that we employed in order to determine the 
value of Kc as described in the previous section. Then 
we determined the value rj[Kc) by doing an extrapolation 
from the effective exponent 77 obtained from the data with 

Lrnax /2 and Lraax- 

The values of 77 that we obtained are reported in Ta- 
ble 1. We also report the value obtained for cr = 1.75 on 
the square lattice and this value extrapolates well with 
the ones on the triangular lattice. In this table, we can 
see that there is near no dependance in the size Lmax- 
This confirms that our determination of Kc is not affected 
by further corrections. 

From Fig. [H we see that in the classical regime and 
in the intermediate regime up to cr ~ 1.5 our results are 
in agreement with the prediction rj = 2 — a (represented 
by a dotted line), confirming the results of Luijten and 
Blote For cr > 2, ry is in perfect agreement with the 
value for a short range model rjsr = 1/4 (represented also 
by a dotted line). In the remaining part for 1.6 < cr < 2 
and shown in the inset of Fig. [1] our results do not agree 
with the prediction of the RG analysis 0, On the 
contrary, we observe a set of values which interpolate 
smoothly between these two behaviors. The errors that 
we obtain on the measured values of 77 are much smaller 
than the deviation observed to the form max{2 — a, rjsr). 
In particular, for cr = 1.8, the deviation corresponds to 
near 10 times the standard deviation, i.e. rj — rjsr — 
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TABLE L Exponent 77 for cr = 1.4-2.2. The second to fourth 
column contain the values obtained by a fit with data up to 
linear size L = 1280, 2560 and 5120. The last column contains 
the predictions from the RG analysis, rj = max{2 — a, 1/4). 



10^77. In the inset of Fig. [l] we also show the results of 
Luijten and Blote from [7|. We note that their results 
are compatible with ours but since our error bars are one 
order of magnitude smaller than theirs, we are able to 
observe a deviation from the RG predictions. In fact, we 
can observe that Luijten and Blote results are compatible 
with both the RG predictions and our results 13 1. 

In this letter we obtained strong evidences that for 
the long range interaction Ising model in 2d, the expo- 
nent 77 does not behave as predicted by RG studies in 
0, 0] , in particular in the intermediate regime and close 
to the boundary with the short range regime as shown 
in Figdiand in Table 1. We believe that further studies 
are needed to recheck the RG analysis. In particular, it 
seems that the wave function normalization is not trivial 
for the long range interaction and this could then give 
further contributions which were not taken in account in 
the previous studies [Hj. 

I thank M. Rajabpour for many valuable discussions 
and advices and T. Blanchard for a critical reading of 
the manuscript. 
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